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Abstract. A new spectral code for constructing general-relativistic models of rapidly rotating stars with an 
unprecedented accuracy is presented. As a first application, we reexamine uniformly rotating homogeneous stars 
and compare our results with those obtained by several previous codes. Moreover, representative relativistic 
examples corresponding to highly flattened rotating bodies are given. 
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1. Introduction 

The study of relativistic, axisymmetric and stationary, 
uniformly rotating perfect fluid bodies is motivated by 
extraordinarily compact astrophysical objects like neu- 
tron stars. Several numerical codes have been devel- 
oped in order to calculate the structure and the gravita- 
tional field of these bodies (Bonazzola & Schneider 
Wilso n |197 2|, But terworth & Ipser 
et al 



1975 



198e, 198C Lattimer et al. 1990 



Herold 1992, Hcrold & Neugebauer 



1989a, 1989b, Eriguchi etal. 1994, Stergioulas & Friedman 



1974 



1976, Friedman 



Neugebauer & 
19921, Komatsu et al. 



1995, Bonazzola et al. 1993; for reviews see Friedman 1998 



and Stergioulas 1998). While they obtain an accuracy of 
up to 5 digits for sufficiently smooth equations of state, 
these methods yield fewer than 4 digits in the case of stiff 
equations of state (e.g., for constant density), which is 
due to particular Gibbs phenomena at the star's surface. 
In order to avoid these Gibbs phenomena, Bonazzola et al. 
(1998) used a multi-domain spectral method with which 
they were able to achieve an accuracy of 12 digits for the 
Maclaurin sequence of homogeneous Newtonian bodies. 

In this Letter we introduce a new numerical code, 
which is based on a multi-domain spectral method for 
representing all metric functions. We intend to use this 
method to investigate neutron stars with realistic equa- 
tions of state. In particular, our multi-domain method 
lends itself to considering several layers inside the star, 
which are characterized by different equations of state. 
As we will outline below, we obtain a hitherto unobtain- 
able accuracy which permits its application even in lim- 
iting cases such as the mass-shedding limit. Moreover, 
we are able to study extremely flattened, homogeneous 
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Einsteinian bodies. Such bodies were the subject of inter- 



esting conjectures made by Bardeen (1971). 

As a flrst application of our method, we reexamine a 
particular example of a uniformly rotating homogeneous 
star that was used by Nozawa et al. (1998) to compare 
three different codes. We give our results, which possess 
a substantially higher accuracy. Moreover, we discuss rep- 
resentative, relativistic examples corresponding to highly 
flattened rotating bodies. 

In what follows, units are used in which the velocity of 
light as well as Newton's constant of gravitation are equal 
to 1. 



2. Metric tensor and field equations 

The line element for an axisymmetric, stationary, uni- 
formly rotating perfect fluid body assumes in Lewis- 
Papapetrou coordinates (p, ^, ip, t) the following form: 



W^dif^] 



e^^{dt + adipy 



To deflne the coordinates (p, C,) uniquely, we require that 
the metric coefficients and their first derivatives be conti- 
nous at the surface of the body. 

In the vacuum region, there emerge three field equa- 
tions of second order to determine the potentials C/, a and 
W . The function k follows from the other potentials by a 
line integral]^. 

In the interior of the body we use the metric functions 
valid in the comoving frame of reference. Here, the only 
new coordinate is Lp' = (p — VLt, where VL is the angular 
velocity of the body. The corresponding line element also 



^ In the highly relativistic regime, we use e^*^, ae'"^ and k — U 
instead oiU ,a and k in order to avoid problems in the presence 
of ergoregions. 
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assumes the above form with potentials C/', a', k' and W , 
which are given by 

e =e [(l + ilaj — it W e J, 
(1 - na')e^^' = (1 + Qa)e^^, 
k' ~U' ^k-U and W ^ W. 

Since in the comoving frame the energy-momentum tensor 
reads 



ik 



4 ; 



where is the mass-energy density and p the pressure, the 
field equations assume a particularly simple form, see, e.g., 
Kramer et al. (1980), equations (19.35a-c). For a given 
equation of state, p = or = ^(p), the relativistic 
Euler equations T^^ -fc = yield 



exp 



dp' 



^^{p')+p' 



= const. 



Hence, pressure and density can be expressed in terms of 
U' . At the surface B of the star, the pressure vanishes 
which leads to a constant surface potential U' — Vq. In 
particular, for homogeneous stars (/i = const.) we obtain 



p = n{e 



Vo-U' _ 



Taking formula (19.35d)0 and the condition (19.37) of 
Kramer et al. (1980) we may express k' by the poten- 
tials a', U' and W via a line integral. Thus again a system 
of three field equations emerges. 

All potentials satisfy regularity conditions at infinity 
and along the axis of rotation (p = 0) and possess more- 
over reflectional symmetry with respect to the plane C = 
(see Meinel & Neugebauer 19951 ). 



3. The numerical scheme 

The numerical scheme to solve the field equations with 
respect to boundary and transition conditions is based 
on a multi-domain spectral method. After imposing re- 
flectional symmetry, the set of all relevant (p, C)~values, 
{(p, C): yO > 0, C > 0}, is divided into several subre- 
gions for physical reasons. In the simplest case, we only 
take two subregions, the interior and the exterior of the 
body. However, if we consider several layers inside the star, 
which are characterized by different equations of state, we 
will be forced to allow for more than two regions. In this 
Letter, we will restrict ourselves to only two subregions. 

Each of these subregions is mapped onto the square 
P = [0, 1]^. In order to do this we introduce a function 
ys defined on the Interval / = [0,1] with 

2/b(0) = 1, z/b(1) = 0, 

^ Note that the bar over C, in equation (19.35d) of Kramer et 
al. (L980) is a misprint. 



which describes the surface of the body by 

S-{(p,C): ^rlt.e ^rlvBit), < i < 1}, 

where and Vp are the equatorial and polar coordinate 
radii of the body respectively. As a particular example 
for the mapping in question, we used for the interior the 
transformation 



P 



and for the exterior 



«2 



In this manner, the axis p — and the plane C = are 
mapped to the coordinate boundaries t — and t — 1 
respectively. Furthermore, the surface B of the body is 
characterized by s = 1. For the interior and exterior trans- 
formation, the point s = corresponds to the origin and 
to infinity respectively. 

We assume all potentials to be smooth functions on 
such that we may approximate them well by two- 
dimensional Chebyshev-expansions with respect to the 
coordinates s and t. In the same manner, we represent 
the unknown boundary function j/s as well as the bound- 
ary values a'g and Wb of the potentials a' and W in 
terms of (one-dimensional) Chebyshev-polynomials with 
respect to the coordinate t. If these three one-dimensional 
functions were given, we would have to solve a particular 
interior and exterior boundary value probleixi^ of the re- 
spective field equations. However, we have to deal with 
a free boundary value problem, where these three func- 
tions are not known, but have to be determined such that 
the normal derivatives of the potentials U', a' and W be- 
have continuously at the surface B ^. Taking only a fi- 
nite number of Chebyshev-coefficients into account for the 
interior potentials U', a' and W, the exterior potentials 
U, a and W and the surface quantities ys, a'g and VFb, 
our numerical scheme consists in determining all unknown 
Chebyshev-coefficients by satisfying the interior and ex- 
terior field equations at a number of grid points in P and, 
moreover, requiring the above transition conditions at the 
surface. The total number of unknown coefficients equals 
the total number of equations. The system is solved by a 
Newton-Raphson method, where the initial guess for the 
entire solution in the case of constant density can be taken 
from the analytical Newtonian Maclaurin solution. 

4. Results 

4.1. Comparison with previous codes 



Nozawa et al. (1998) compared three different codes for 
various choices of the equation of state. We take one par- 
ticular example in which for constant density they pre- 
scribed the normalized central pressure Pc — Pc/p — 1 

^ The boundary values are completed by U' = Vq. 
* It is a consequence of the field equations that k' is then 
also difi'erentiable. 
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(1) 


(2) 


(3) 


Pc 


1 








rp/re 


0.7 






0.11% 


n 


1.41170848318 


1.1% 


0.32% 


0.97% 


M 


0.135798178809 


2.3% 


0.19% 


0.86% 


Mo 


0.186338658186 


0.17% 


0.32% 


1.4% 


Rcirc 


0.345476187602 


0.098% 


0.053% 


0.27% 


J 


0.0140585992949 


1.6% 


0.045% 


2.3% 


Zp 


1.70735395213 


6.1% 


0.013% 


2.1% 




-0.162534082217 


1.7% 


1.9% 


4.4% 


^eq 


11.3539142587 


17% 


0.10% 


8.1% 



Table 1. Results for the constant density model with pc — 
1, fp/re = 7 X 10^^, and relative errors of the codes (1): 
KEH(OR), (2): KEH(SF) and (3): BGSM. ^ = n/n^'"^, 
M = Mfi^/^, Mo = MoA^^/^ Rdrc = RdrcH^^'^ and 
J = J ^ are normalized values of the angular velocity SI, 
gravitational mass M, rest mass A^o, equatorial circum- 
ferential radius Rdrc and angular momentum J. Zp is the 
polar redshift, Zl^ and Z^^ are equatorial redshifts for 
photons emitted in the forward and backward direction. 



and the ratio rp/re = 7 x 10""'^. The results of this com- 
parison are given in table 11 of Nozawa et al. ( 1998| ). 
Here we calculate the same quantities with an accu- 
racy of 12 digits and list them in column 2 of Table 
Columns 3-5 refer to the codes by Komatsu et al. 
1989h| ) and Eriguchi et al. ( [1994D [abb reviated 



ia,c 



(|l989a 



by KEH(OR)], Stergioulas and Friedman (|1995| ) [abbre- 



viated by KEH(SF)] and Bonazzola et al. ( |1993D [abbrevi 
ated by BGSM] and give the relative error of the quantity 
in question, i.e. for example, |M[KEH(OR)] - M\/M « 
0.023. The fact that Vp/re was not exactly 0.7 in the 
BGSM calculation does no t affe ct the comparison sub- 
stantially. In Nozawa et al. ( 1998 ), the general-relativistic 
virial identitie s GR V2 and GRV3 (derived by Bonaz zola fc 
Gourgoulhon 1994 and Gourgoulhon & Bonazzola 1994 ) 
were calculated to check the accuracy of the numerical 
solution. For our result, this check yields 1.6 x 10~^^ for 
GRV2 and 4.4 x 10"" for GRV3. Note that we used 23 
Chebyshev-polynomials for each dimension in this calcu- 
lation. As an additional test of accuracy, we calculated 
the angular momentum and the gravitational mass in two 
different ways: (i) from the asymptotic behaviour of the 
metric and (ii) by means of integrals over the matter distri- 
bution [cf. Bardeen & Wagoner ( 1971 ), equations (11.24), 
(11.26) and (11.23), (11.25) respectively]. We get a relative 
deviation of 1.7 x 10^^"' for the mass and 6.2 x 10"^'' for 
the angular momentum. 



^ Note that there are misprints in Nozawa et. al. (1998) con- 
cerning the formulae for Zp, Zlq and Z^^: In equation (26), 
—2up has to be replaced by —Vp, and in equations (27) and 

(28), (i^e -/3e)/2 by (/3e -i^e). 





(a) 




(b) 




(c) 




Pc 


0.002 




0.004 




0.003 




Tp/re 


0.2 




0.1 




0.1 




Q. 


1.089468 e+00 


1.004798 e-fOO 


9.305304 e- 


-01 


M 


8.314919 e- 


-04 


1.099261 e- 


02 


5.603298 e-^ 


-03 


Mo 


8.353271 e- 


-04 


1.124879e- 


-02 


5.677891 e- 


-03 


Rcirc 


1.028186 e- 


-01 


2.474602 e- 


-01 


2.249448 e- 


-01 


J 


3.663135e- 


-06 


3.093789 e- 


-04 


1.072577e- 


-04 


Zp 


1.585397 e- 


-02 


9.132868 e- 


-02 


5.381392 e- 


-02 


Ztq 


-9.890231e- 


-02 


-1.912324 e- 


01 


-1.720749 e- 


-01 


yb 

^eq 


1.308359 e- 


-01 


3.827610e- 


-01 


2.825350e- 


-01 


R 


0.2444465 




0.3001205 




0.3522896 




GRV2 


l.le-09 




3.1e-ll 




1.9e-08 




GRV3 


1.6e-09 




3.2e-ll 




2.3e-08 





Table 2. Three highly flattened, constant density models. 
The corresponding surface shapes can be found in Fig. 1. 
In these calculations we used 23 Chebyshev-polynomials 
for each dimension. 



4.2. Highly flattened rotating bodies 



Chandrasekhar (1967) has shown that the post- 
Newtonian Maclaurin spheroids become singular at the 
eccentricity e — ei — 0.98522. . . , the poi nt of t he first ax- 
isymmetric secular instability. Bardeen ( |l97l|) confirmed 
this result and discussed the possibility of two Newtonian 
axisymmetric sequences bifurcating from the Maclaurin 
spheroid at e = ei, a first one that should finally evolve 



son 


189^ 




L89J 


), Wong 


1885a 




L885b 


1885c 



see also Lichtenstein 1933) and a second one, which he 
called the 'central bulge configuration', that sh ould end 
in a mass-shedding limit. Eriguchi & Sugimoto ( |l98lD in- 
deed found the first sequence, and called it the 'one-ring 
sequence'. The properties of the bifurcation have been an- 
alyzed by Christodoulou et al. ( |1995D . 

By means of a Newtonian version of our code we 
were able to find the second sequence as well. It ap- 
pears as a continuation of the one-ring sequence across the 
Maclaurin sequence and is characterized by a (bi-convex) 
'lens shape' of the solutions. 

The point e = ei of the Maclaurin sequence corre- 
sponds to the value R = Ri = 0.27320. . . of the rota- 
tion parameter R = J'^ /Mq^^^. Bardeen ( 1971 ) speculated 
that there should be a gap in the i?-values of relativis- 
tic solutions around i?i and that the relativistic solutions 
should show some properties of the Newtonian 'central 
bulge' or Dyson ring solutions when approaching this gap 
from the sphere end or from the disc end of the Maclaurin 
sequence. 

In order to check these conjectures, we have calculated 
three relativistic models in the vicinity R — Ri. These 
models are characterized by parameters as shown in Table 
2 and their (coordinate) shapes are depicted in Fig. 1. Note 
that the high accuracy of our code is crucial when investi- 
gating the subtle behaviour of the relativistic solutions in 
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(b) 












p 







Fig. 1. Meridional cross-sections of the solutions specified 
in Table 2. 



this region — an earlier attempt by Butterworth ( 1979 ) 
could not finally clarify these questions. Solution (a) shows 
indeed a 'lens shape', whereas solution (b) has a 'one-ring 
tendency' as does solution (c), albeit far less pronounced. 
The solution (a) is rather close to the Newtonian 'lens 
shape' sequence. The situation becomes more and more 
complex when the bifurcation points R2 = 0.36633. . . , 
i?3 = 0.43527. . . , etc. (corresponding to €2 = 0.99375. . . , 
€3 = 0.99657. . . , etc.) of the Newtonian two-ring, three- 
ring, etc. sequences are approac hed (f or the two-ring se- 
quence see Eriguchi & Hachisu 1982| ). Our solution (c) 
already has a value for R close to i?2, see Table 2. 

A detailed analysis of highly flattened, rotating 
Newtonian as well as Einsteinian bodies including a 
discussion of stability aspects and of the route to in- 
finitesimall y thin, relativistic disc s (Ba r deen & Wagoner 



1969| , |1971| , Neugebauer & Meinel |1993| , |1995| ) will be the 



subject of a future publication. 
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